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Abstract 

In a freely cooling granular material fluctuations in density and temperature cause position 
dependent energy loss. Due to strong local dissipation, pressure and energy drop rapidly 
and material moves from 'hot' to 'cold' regions, leading to even stronger dissipation and 
thus causing the density instability. The assumption of 'molecular chaos' is valid only in the 
homogeneous cooling regime. As soon as the density instability occurs, the impact parameter 
is not longer uniformly distributed. The pair-correlation and the structure functions show 
that the molecular chaos assumption — together with reasonable excluded volume modeling 
— is important for short distances and irrelevant on large length scales. 

In this study, the probability distribution of the collision frequency is examined for 
pipe flow and for freely cooling granular materials as well. Uncorrelated events lead to a 
Poisson distribution for the collision frequencies. In contrast, the flngerprint of the cooper- 
ative phenomena discussed here is a power-law decay of the probability for many collisions 
per unit time. 

Keywords: discrete element method, event driven simulations, clustering instability, arch- 
ing, shock waves, power-law distribution, cooperative phenomena. 
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1 Introduction 

Many rather astonishing phenomena are known to occur when granular materials like sand 
or powders move |l|-@]. Of interest are, e.g. density waves emitted from outlets P], crack 
formation during vibration or during the flow through a pipe l^j-Q, and pattern formation 
due to dissipation []3], piO]-]!^ . AH these effects are connected to the ability of granular 
materials to form a hybrid state between a fluid and a solid: energy input can lead to a 
reduction of density so that the material becomes 'fluid' and, on the other hand, in the 
absence of energy input, granular materials 'solidify' due to dissipation. Thus, a packing of 
sand behaves like a solid when pushed, but offers no resistance to a pulling force. 

In order to formalize and quantify the complicated rheology of granular media various 
attempts have been made. Continuum equations of motion and kinetic theories p^]-|20[ 



are 



the first successful steps towards a quantitative description of granular materials — at least 
for limited parameter range. This restriction exists because it is very difficult to incorporate 
into these theories all the details of static friction or other relevant microscopic mechanisms. 



Also the generalization to high densities is an ardous task, see for example Refs. pT],p2 
and references therein. Most of the classical, but also the more advanced theories are based 
on the assumption of molecular chaos — the assumption that the velocities and the relative 
positions of all colliding pairs of particles in a gas are uncorrelated. In a dilute gas, the errors 
introduced by this assumption are negligible. In dense granular ffows, however, correlations 
between colliding particles may be important, leading to qualitative changes of behavior. 

Section || is dedicated to a brief introduction of the different discrete modeling ap- 
proaches. In particular, we will present the 'hard-particle' Event Driven (ED) [p3],p^ and 
the 'soft-particle' Molecular Dynamics (MD) |2^, 2^, 26] methods, both for inelastic spher- 
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ical particles with frictional forces. We extend the traditional ED method by introducing 
a cut-off time tc for dissipation |1^,|2^,^. Any particle that encounters a second collision 
before this time passed by is assumed to be elastic — the extended method is named TCED 
| |28| . Furthermore, the Direct Simulation Monte Carlo (DSMC) approach is discussed and 
applied to freely cooling granular media. The validity of the molecular chaos assumption 
in granular flows was examined by comparing event driven (ED) 'hard sphere' simulations 
to those performed with the Direct Simulation Monte Carlo (DSMC) method [Q]. The ED 
method is capable of reproducing velocity correlations — even in the limit of rather large 
densities — whereas DSMC assumes molecular chaos by construction. In Section |^ the 
structure factor and the pair-correlation function are examined and reasons for the break- 
down of molecular chaos are discussed. The clustering instability is described in section ^ 
with respect to the restitution coefficient and a cut-off time for dissipation. The probability 
distribution function for the collision frequency is measured in the homogeneous and the 
non- homogeneous clustering case and differences are evidenced 0|. The same is also found 
in pipe flow ^ where shock waves and arching are the observed cooperative phenomena. 
The probability distribution functions are measured and discussed in Section | and interest- 
ingly have the same functional behavior as in the case of the clustering instability. Finally, 
the results are summarized in Section 0. 



2 Models for particle-particle interactions 

The basic constitutents of granular materials are mesoscopic grains, made of, for example, 
10^° molecules. When these objects interact (collide) the attractive potentials of the indi- 
vidual atoms can often be neglected. Three models for the particle-particle interactions are 
discussed in the following. They account for the excluded volume of the particles via a re- 
pulsive potential, either 'hard' or 'soft', or assume point particles and introduce appropriate 
corrections. 

It is important that the surface of the grains is rough on a microscopic scale so that 
solid friction occurs. In general, one has to distinguish between sliding, sticking, and rolling 
friction, but we will only discuss simplified models here. An entire discipline called tribology 
has evolved to study solid friction in depth [|l], ^ . Friction and other sources of dissipation, 
like viscous damping or plastic deformations, have the crucial consequence that the system 
does not conserve energy. Since dissipation may occur due to various reasons, we discuss 
in the following only simple dissipation laws, assuming that the detailed knowledge of the 
interaction potential is of minor importance. In fact, complicated laws often increase the 
number of parameters without giving qualitatively different results |^ . 



The difference between the two most frequently used discrete element methods is the 
repulsive interaction potential. For the molecular dynamics (MD) method, soft particles 
with a power-law interaction potential are assumed, whereas for the event driven (ED) 
method perfectly rigid particles are used. The consequence is that the duration of the 
contact of two particles, tc, is finite for MD, but vanishes for ED. In the DSMC method, one 
assumes point particles without repulsive potential but with an effective scattering cross 
section. In addition one applies corrections from the kinetic theory, in order to account for 
the effective free volume, the reduced mean free path, the increased collision frequency, and 
the enhanced collisional momentum transport. 
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2.1 The event driven, rigid particle method 

Here, we apply the simplified collision model introduced by Walton and Braun and 



recently experimentally established by Foerster et al. |]32[ and Labous and Rosato |3^. For 
given velocities before contact, three coefficients are needed to evaluate the velocities after 
collision. At first, the coefficient of normal restitution, r, defines the incomplete restitution 
of the normal component of the relative velocity. The second, the coefficient of friction, /i, 
relates the tangential momentum change to the normal one, i.e. Coulomb's law. The third, 
the coefficient of maximum tangential restitution, /3o, delimits the restitution of tangential 
velocity of the contact point to ensure energy conservation. Note that this model implies that 
two grains at contact either slide, following Coulomb's law, or stick together [pT|,p2tp^, ^ 



In the following, we apply the basic conservation laws and determine the equations for the 
velocities after collision. 

Consider two particles with diameter di and d2 and masses mi and m2. The normal 
unit vector for their contact is n = {fi — / — 'r2|, where rj is the vector that gives the 
position of the center of particle i (i = 1, 2). For the interaction of particle i = 1 with a 
fixed wall, we set m2 = oo, and n is in this case the unit vector perpendicular to the wall 
surface pointing from the contact point with the wall to the center of the particle. The 
relative velocity of the contact points is 

Vc = vi-V2- + ^ ^ ' (1) 

with Vi and Coi being the linear and angular velocities of particle i just before collision. From 
the momentum conservation laws for linear and angular momentum follows 

v'^ = vi + AP/mi (2) 

^[=^1- j^n X AP , (3) 

v'2 = V2 — AP/m2 , and (4) 

^'2 = ^2- j^n X AP , (5) 

where and a;,' are the unknown velocities of particle i after collision. Jj is the moment of 
inertia about the center of particle i and AP is the change of linear momentum of particle 
1 and is a function of r, fi, and Pq: 

AP = -mi2(l + r)t;(") - ^mi2(l + f3)vi'\ (6) 

with the reduced mass mi2 = 77111112/ {nii + 7712). {71) and (t) indicate the normal and the 
tangential components of Vc, respectively, and the factor 2/7 in the tangential part of Eq. (|^) 



stems from the fact that solid spheres are used |^ . r is the coefficient of normal restitution 
and /3 = min[/3o,/3i] is the coefficient of tangential restitution. The latter is simplified so 
that either sticking or sliding are exclusively allowed. A sticking contact has a constant 
maximum tangential restitution f] = Pq, with — 1 < /?o ^ 1 due to the elasticity of the 
material. Typical values for, e.g. acetate or glass, are /?o ~ 0.5 Sliding, Coulomb-type 
interactions have /3 = i.e. AP^*) is limited by jU^AP*^"^. Using the basic conservation 
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laws one can find (3i = —1 — + r) cot(7)(l + l/qi), with the colhsion angle 7, and the 
factor Qi = 4Jj/ {rriidf) that accounts for the mass distribution inside the particles [p2|p^, p5[ . 
As illustration, a schematic picture of two colliding particles is given in Fig. [l|. The angular 
velocities are ui = 002 = immediately before collision (a) and non-zero after collision (b). 
For a more detailed discussion of the above equations see Ref. [0, ^ . 



For the simulation of rigid particles, we use an event driven method such that the 
particles undergo an undisturbed motion until an event occurs. An event is either the 
collision of two particles or the collision of one particle with a wall. From the velocities 
just before contact, the particle velocities after a contact are computed following Eqs. (0)- 
d^). Lubachevsky introduced an efficient scalar ED algorithm which updates only those 
particles involved in the previous collision. Like in Refs. ^] we implement the algorithm 
of Ref. p4[ with some changes and extensions. Despite gravitational acceleration, all times 
of contact of particles with each other or with the lateral walls can be calculated analytically. 
The coefficient of normal restitution depends on the partner of the colliding particle, i.e. r 
or is used to indicate particle-particle, or particle-wall collisions, respectively. 



2.2 The connection between hard- and soft-sphere models 

In the ED method, the time during which two particles are in contact tc is implicitly zero. 
The consequence is that exclusively pair contacts occur and the instantaneous momentum 
change AP in Eqs. (0)-® suffices to describe the collision completely. ED algorithms with 
constant r run into difficulties when the time between events, t^v, becomes too small — 
typically in systems with strong dissipation — and the so-called 'inelastic collapse' occurs 
36|- p^] . To handle this problem, several attempts have been proposed recently [Il2|,p7|p8 



5^]-^^, and one of them |]T^^,^| switches off dissipation when the collision frequency 



becomes too large. 

In MD simulations of dynamical systems, on the other hand, tc > and only a limited 
amount of energy can be dissipated per contact duration. A finite tc thus implies a finite 
energy dissipation rate. In a dense system of soft particles, energy dissipation becomes 
ineffective, i.e. the 'detachment effect' occurs |^,^. This effect is not obtained with hard 



particles and a constant coefficient of restitution r, however, the effect can be observed when 
using 



r for t« > tc . . 

1 for tW<t, , 



as the restitution coefficient for the collision n of particle i. In Eq. (^) t^^ is the time 
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since the last collision and tc is the threshold for elastic contacts that can be identified (up 
to a constant factor of order unity) with the contact duration in the soft particle model. 
Thus, an additional material parameter is defined for the hard sphere model, that leads to 
qualitative agreement between ED and MD simulations and in addition avoids the inelastic 
collapse. Note that tc has different physical meaning in either hard or soft sphere model. 
The traditional ED method has tc = 0. The extended model that uses Eq. (|^ with tc > 
is in the following referred to as the TCED method. 

The integral of all forces f(t), acting on a particle at time t G [to, tc], is needed to 
calculate the momentum change of this particle in the framework of the soft particle model: 

rto+tc _, 

AP = f{t)dt. (8) 

In general, the contact begins at time to and ends at time to + tc- For a constant force / or 
an infinitesimally small time interval tc, the momentum change AP in Eqs. (@)-(0) can be 
replaced by the term f{t)dt to arrive at a differential formulation for the change of velocities 
v' — V and uj' — uj. The primed and unprimed quantities are the values at time to + ic and 
to, respectively. 

The stress tensor, defined for a test volume V^, can be written as 



(9) 



The indices a and /3 are the Cartesian coordinates, the Tq, are the components of the vector 
from the center of mass of a particle to the point j, where a force with components fp acts. 
Particle i has the mass rrii and a velocity with the components Va- The first sum runs over 
all contact points j, and the second sum runs over all particles i, both within Vc- In the 
static limit, the second term vanishes, since all velocities vanish. On the other hand, for a 
hard-sphere gas, the left term has to be treated differently, since no forces are defined. The 
dynamic equivalent to fp is the change of momentum per unit time AP^/At. For a hard- 
sphere gas the stress due to collisions may be evaluated as the average over all collisions in 
the time interval At. The first sum runs over all collisions taking place in the time between 
t — At and t. In general, the volume Vc and the time-interval At have to be chosen large 
enough to allow averages over enough particles and enough collisions. 



2.3 The time driven, soft particle technique 

Even without using the soft particle method in this study, it is convenient to discuss briefiy 
the standard interaction forces and their connection to the hard sphere collision operator 
that involves the total momentum change AP. Replacing AP in Eqs. (0)-(^ by /(t)tMD, 
with the molecular dynamics time step tMD, allows the integration of the corresponding 
equations of motion with standard numerical methods p3|,p6|. 

Since the modeling of realistic deformations of the particles would be much too com- 
plicated, let us assume that the overlap of two particles is the quantity important for the 
interaction potential. The interaction is short range, i.e. the particles interact only when 
they are in contact so that their penetration depth 5 = \{di + c?2) — (ri — ■ n is positive. 

The first force, acting from particle 2 on particle 1 — accounting for the excluded 
volume which each particle occupies — is an elastic repulsive force 



(10) 
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where /c„ is the elastic modulus and Sq is a normalization constant dependent on the non- 
linearity u and the dimension. In the simplest case of a linear spring that follows Hooke's 
law, u = 1, in the case of elastic spheres in three dimensions, u = 3/2, i.e. a Hertz contact 
and for conical contacts, z/ = 2 can be used. 

The second force — accounting for dissipation in the normal direction — is a viscous 
damping force 

/diss = lJi6/6o)'^n , (11) 

where 7„ is a phenomenological viscous dissipation coefficient and 6 = —vu-n = —{vi—V2)-n 
is the relative velocity in the normal direction. 

The simple linear spring-dashpot model (with u = 1 and = 0) can be solved analyti- 
cally and leads to a contact duration tc = n/u and a restitution coefficient r = exp(— Trrz/u;), 



with u = JuJq - rf, ul = kn/mu, r] = 7„/(2mi2), and rriu = mim2/(mi + 1712) A 



nonlinear repulsive potential can at least be solved in the limit 7„ — and the dependency 



of r on the velocity of impact can be estimated with reasonable accuracy ||46|,^,|50 . 

The third force — accounting for friction in a simplified way — acts in the tangential 
direction and can be chosen in the simplest case as 

/shear = "7*^^ , (12) 

where jt is the viscous damping coefficient in tangential direction and C, = vu ■ t is the 



tangential component of the relative velocity, with t = Vu/lvul- Eq. (JV4 ) is a rather 
simplistic description of shear friction. For many applications (arching, heap formation) it 
is, however, important to include more realistic static friction ||l],0 what can be realized 



by a virtual tangential spring p5|,|5^: when two particles start to touch each other, one 
puts a virtual spring between the contact points of the two particles, and ^{h) = Jl^^[vi2{t) ■ 
i{t)]dt i(ti) is the total tangential displacement of this spring during the contact. The 
restoring frictional force is thus —kf^. According to Coulomb's criterion, the maximum 
value of the restoring force is proportional to the normal force /^ at this contact, with the 
friction coefficient /i. Cast into a formula this gives a friction force 

^ric = -4min(A;,|e|,/i/;;) . (13) 

Note that the tangential spring has to be kept at a maximum length ^max = f^fn/^t in 
order to lead to reasonable agreement with contact dynamics simulations or theoretical 
calculations [Q. Only when particles are no longer in contact with each other is the spring 
removed. The main source of static friction in real systems is the geometrical roughness of 
the surfaces [p4|-^, and the same effects of particle stopping can be obtained also without 



Eq. (|T3D by using particles of complicated shapes, like crosses or polygons ||58|-|6T|. In fact 



when particles deviate from the spherical shape, rotations are suppressed in dense packings 
under strong load. However, in some cases it is sufficient to use a combination of Eqs. (0) 
and dll): 

^yn = - min(7j^, fif^)t , (14) 
a rather bold alternative to the more realistic static friction law in Eq. (|T3D, but sufficient 



for many, especially dynamic situations [Q. For a comparison of ED and MD simulation 
results, see Refs. P,ra,|5l,lK 
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2.4 DSMC simulation method 

Direct simulation Monte Carlo (DSMC) was first proposed for the simulation of rarefied 
gas flows [0 and is also used for liquid-solid flow simulations, see Ref. and references 
therein. One of the algorithm's advantages is its suitability for parallelization, what makes 
it a convenient tool for the modeling of granular media p5|,|65|. 

In DSMC the evolution of the system is integrated in time steps r, at each of which 
every particle is first moved without feeling other particles. The particles are then sorted 
into square cells with sides Lc and volume Vc = L^. L^. is set as one half of the mean 
free path, but never smaller than two bead diameters. The time step r is always chosen 
small enough to assure that even the fastest particle needs several time steps to cross a 
cell. Between particles in the same cell stochastic collisions take place; the rules for these 
collisions are taken from kinetic theory. First we choose the number of collision pairs in 
each cell, 

Mc = — / , (15) 

where A^^ is the number of particles in the cell, Vmax is an upper limit for the relative velocity 
between the particles, and a = 2d is the scattering cross section of discs in 2D. To get v^ax 
we sample the velocity distribution from time to time and set v^ax to be twice the maximum 
particle velocity found. In order to account for the actual relative velocities we apply an 
acceptance-rejection method: for a pair of particles i and j the collision is performed if 

\Vi - Vj\ 



< Z , (16) 

Umax 

where Z is uniformly distributed in the interval [0.0, 1.0[. This method leads to a collision 
probability proportional to the relative velocity of the particles. 

Since the collision takes place regardless of the position in the cell, we have to choose an 
impact parameter h in order to calculate the post-collision velocities. The impact parameter 
is defined as 



\Vi - Vj\ 



(isin7 , (17) 



where 7 is the angle between [vi — vj) and (fi — fj). For central collisions 6 = 0, and b = d for 
grazing collisions. Assuming molecular chaos, b is drawn from a uniform distribution in the 
interval [— c?, c?]. The rest of the collision scheme is identical to the event driven procedure, 

so that the normal component of the post-collision velocity is v'^ ^ = —rv^^\ whereas the 
tangential component remains unchanged, i.e. jSo = —1. 

To get better results at higher densities the DSMC method was modified in two re- 
spects. First, the number of collisions Mc in equation (|15D is increased by replacing the 
volume Vc of a cell with the effective free volume — Vq, where Vq is the volume the par- 
ticles in that cell would need in a random close packing with packing fraction 0.82 in 2D 
66| . Second, we added an offset of d to the particle distance along the direction of the 



momentum transfer after the collision [H^ . 



3 Comparison of ED and DSMC simulations 

In this section two simulations are presented, starting with the same initial condition, and 
using the same parameters, but being carried out with the ED and the DSMC methods. 



3 COMPARISON OF ED AND DSMC SIMULATIONS 



9 




Figure 2: Normalized kinetic energy vs. normalized time from an ED and a DSMC 
simulation in 2D with identical initial conditions and N = 99856, g = 0.25, r = 0.8. The 
dashed line represents Eq. (|TBp. The 

In ED the probability distribution of the impact parameter may deviate from the case 
expected for molecular chaos, whereas DSMC always uses a constant probability for b in 2D. 
The simulations involve N = 99856 = 316^ dissipative particles with restitution coefficient 
r = 0.8 in a periodic quadratic system with volume fraction g = 0.25. The system size is 
I = Ld with dimensionless size L and particle diameter d. In order to reach an equilibrated 
initial condition, the system is first allowed to evolve with r = 1 for about 10 collisions 
per particle, so that a Maxwellian velocity distribution and a rather homogeneous density 
distribution exists. Then, at t = Os, dissipation is set to r = 0.8 and the quantities of 
interest are calculated as functions of time. 



3.1 Freely cooling granular materials 



In the homogeneous cooling state |T0|, |T6|, ^ we expect that the energy K{t) of the 

system decays with time and follows the functional form 



1 



with the theoretically expected time scale 

y/TTds^{Q) 



2{l -r'^)QV 



(19) 



the initial kinetic energy Kq = 
and the volume fraction g with 



0.36, 



as a function of the initial mean velocity v = ^J2Kq/ Nm 
fr(0), the particle diameter d, the restitution coefficient r 

s.{q) = (1-^)7(1-7^/16) nmn 

s^{q) ^ 0.63158, d = 0.001m, and v 

In Fig. we present the normalized kinetic energy K{t)/Ko as a function of the 
normalized time t/to. At the beginning of the simulation we observe a perfect agreement 



Inserting the corresponding parameters 1- 
0.2047m/s, we have Iq^ = 23.24 s'^ 



3 COMPARISON OF ED AND DSMC SIMULATIONS 



10 



5 




Figure 3: Normalized probability distribution of the contact parameter from an ED sim- 
ulation in 2D with N = 99856, g = 0.25, and r = 0.8 at different times, as given in the 
inset. 

between the theory for homogeneous cooling and both simulations. At t/to ~ 2 substantial 
deviations from the homogeneous cooling behavior become evident, and only at t/to ^ 
10 a difference between ED and DSMC can be observed. After that time, the kinetic 
energy obtained from the DSMC simulation is systematically smaller than K{t) from the 
ED simulation. We relate this to the fact that the molecular chaos assumption of a constant 
probability distribution of the impact parameter b is no longer valid. Since dissipation acts 
only at the normal component of the relative velocity, DSMC dissipates more energy than 
ED as soon as the number of central collisions is overestimated . To verify this assumption 
we take a closer look at the impact parameter and its probability distribution in the next 
subsection. 



3.2 The impact parameter 

One basic assumption connected to molecular chaos is a uniform probability distribution 
of the impact parameter. We define P(h/d) to be the probability distribution of h and 
normalize it such that /q^ (1xP{x) = 1. From ED simulations with elastic particles the 
normalized probability distributions P{b/d) = 1 in 2D, and P{b/d) = 2b/ d in 3D, are found 
as expected for the case of molecular chaos. 

The ED simulation of Fig. ^ leads to P{b/d) = 1 for short times only. For larger times 
we observe an increasing (decreasing) probability of grazing (central) collisions. In Fig. ^, 
data of the probability distribution are presented at different times during the simulation 
of Fig. 0. As it is obvious from the data, more and more grazing collisions occur with 
increasing simulation time. Evidently, the assumption of a constant probability distribution 
of the impact parameter is violated. 

One can imagine at least two reasons for the deviation of P{b/d) from the constant 
value. The first is that P{b/d) might be a function of the density, and that due to density 
fluctuations, the form of P{b/d) changes. Thus we calculate P{b/d) in smaller systems with 
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= 240, r = 1, and different volume fractions, ranging from very dilute to extremely 
dense systems. P{b/d) is not sensitive to the density, as long as the collisions are elastic 
35|. Another reason for P{b/d) to deviate from unity might be dissipation. In Fig. ^ the 
restitution coefficient is varied for fixed g = 0.7495. For weak dissipation, i.e. r > 0.9, 
the distribution is homogeneous. For stronger dissipation r = 0.80 we find an increasing 
probability of grazing contacts. 



5 
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Figure 4: Normalized probability distribution P{b/d) from ED simulations in 2D with N = 
240, g = 0.7495 and different r. 

The assumption P{b/d) = 1 is true in elastic systems for arbitrary density. For 
inelastic systems, P{b/d) is constant for sufficiently weak dissipation but depends on b/d 
for strong dissipation. The breakdown of molecular chaos is not due to high density, and 
also dissipation is not the only reason for it, since the dissipation must be strong enough to 
cause the inhomogeneous distribution. The remaining question is: why do we observe this 
increasing probability of grazing contacts? 

Looking in more detail at the simulations in Fig. ^, we observe that the inhomogeneous 
distribution for r = 0.8 is connected to shear motion of the particles, whereas no visible shear 
motion occurs for r > 0.9. The shear motion in the left panel of Fig. |^ can be understood as 
the geometrical reason for the higher probability for grazing contacts, i.e. layers of particles 
shear against each other and therefore touch preferentially with large impact parameter b. 
The two eddies (top left and bottom right) in the left panel are rotating in the same sense 
and lead to small stresses in the interior. Outside, where the velocities are less correlated 
one observes larger stresses and the stress tensor typically has strongly different eigenvalues. 
The more disordered situations in the middle and right panel, in contrast, are connected to 
rather random stresses. 



3.3 The structure factor 

One difference between ED and DSMC simulations is the handling of excluded volume by 
the two methods. While ED handles hard spheres with a well defined excluded volume, the 
DSMC method models point particles and excluded volume is introduced by the approx- 
imations described in subsection As expected, we obtain dramatic differences in the 
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Figure 5: Snapshots from ED simulations in 2D with = 240, g = 0.7495 and different 
r = 0.80 (left), r = 0.90 (middle), r = 0.95 (right) after t = 1 s. The lines give the velocity 
vector scaled by 1, 3, and 7 in the left, middle and right panel, respectively. In the panels 
below the stress is plotted in the principal axis representation following Eq. (^. The average 
is taken in the time interval 0.5 s <t< 1.0 s over all collisions in each cell. The maximum 
eigenvalues o"max, measured in the left, middle, and right panel are 0.894, 0.568, and 1.56, 
respectively. They are measured in arbitrary units and the figures are scaled so that (Jmax 
has the same length. 
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particle-particle correlation function g{r/d) in Fig. where r/d is the particle separation 
in units of particle diameters d. Since ED (solid lines) models hard spheres, g{r/d) = 
when < r /d < 1. For short times one observes g as for an elastic hard sphere gas and 
at larger times g{r/d) shows a rich structure with peaks at 1, 2, 3, ... and multiples of 
\/3/2 (factor > 2), indicating a rather close triangular packing of monodisperse spheres. In 
contrast, the DSMC simulations (dashed lines) show no short range correlations between 
particle positions throughout the whole simulation. 
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Figure 6: Correlation function g{r/d) as obtained from the ED (solid line) and DSMC 
(dashed line) simulations of Fig. ^ The different curves are shifted vertically in order to 
avoid overlap. The solid line at < r/d < 1 gives the baseline g{r/d) = 0. 



The next question is, whether this difference has consequences at greater length scales. 
The formation and growth of large clusters [|1^,|6^,|6^ is quantified by the large wavelength 
modes of g (r/d), or equivalently, by the structure factor S{k) at small wave-number k. We 
calculate S{k) by a direct FFT (fast Fourier transform) of the two-dimensional density. 
Before we apply the FFT we map the particles onto a M x M lattice, where M is the closest 
power of 2 that gives a lattice box size of about one diameter. 

The structure factors obtained by ED are presented in Fig. 0(a) and those obtained 
by DSMC in Fig. ^(b). Different symbols correspond to different times. We observe an 
increase of S{k) for short wavenumbers k < 25, until the structure factor ceases to change 
for t > 20 s. The structure factor agrees reasonably well for both simulation methods, and 
for large enough times it does not change further. This proves that the DSMC simulation is 
capable to reproduce the more realistic but computationally more expensive ED results that 
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Figure 7: (a) Structure factor obtained from the ED simulations of Fig. ^ as function of the 
wavenumber k = l/X, with wavelength A and system size /. (b) Structure factor obtained 
from the corresponding DSMC simulation. 

account for the excluded volume by construction. Even without short-range correlations, 
the information about large wavelengths is well reproduced by DSMC simulations. 



4 The clustering instability 

The essential difference between a granular and a classical gas is dissipation. The resulting 
clustering instability was examined in ID p6|-p9|,[7^-|7i| and in 2D |T0|, |69|, |T5|-[7^] . Cluster 



69 led to 



growth could be described theoretically only in the case of irreversible aggregation [ 
the more general case of reversible aggregation is still an open issue. 

Detailed examination of the inelastic collapse by McNamara and Young 
the picture of different 'phases'. In a periodic system without external forcing exists a 
critical dissipation — connected to volume fraction and restitution coefficient — above which 
clustering occurs and below which the system stays in molecular chaos. In the transient 
regime the system behavior seems to depend on the system size and shearing modes are 
frequently observed. 



4.1 Parameter studies 

In the following we examine periodic systems of length L = l/d m 2D, with the particle 
diameter d, using the TCED simulation method as described in subsection |2.2| . The volume 
fraction is ^ = Nn{d/2iy = {n/4)N/L'^ and the particles are initially arranged on a square 
lattice, homogeneously distributed in the system. First the system is equilibrated with 
r = 1, then the dissipation is switched to the desired value and the simulation starts at 
t = s. Here we use a rather small system with = 784, L = 50, and g ~ 0.25. 

The restitution coefficient and the cut-off time are varied (0.99 < r < 0.20 and 



4 THE CLUSTERING INSTABILITY 



15 



10"^'' s < tc < 10~^s) and the simulation is performed for each parameter set until every 
particle carried out C/N = 1000 collisions. In Fig. ^(a) K{t) is plotted against the mean 
number of collisions per particle C/N for simulations with = 10~®s and variable r. For 
large r and small C/N the energy behaves as K{t) oc exp{—C/N) as can be derived from 
Eq. (|l^) by integration over the product of mean velocity v = K/m and inverse mean 
free path A oc 1/q. With decreasing r the initial slope is larger because more energy is 
dissipated per collision. At larger times, energy decays much slower and K{t) deviates from 
the straight line already for r < 0.95. In Fig. H(b) {K^ — Kz)/ K{t) is plotted against C /N . 
In a homogeneous system without clustering, the value of {K^ — K^) / K{t) fluctuates around 
zero, while values close to unity indicate the 'shearing- mode', i.e. most of the kinetic energy 
can be found in one direction [^, The deviations from the homogeneous cooling state 
begin earlier with decreasing r. 
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Figure 8: (a) K{t) as function of C/N with tc = 10 ^s and r as given in the inset, (b) 
{Kx — Kz)/K{t) for the same simulations as in (a). 




Figure 9: (a) K{t) as function of t with r = 0.6 and — log^gtc as given in the inset, (b) 
K{t) plotted against C/N for the same simulations as in (a). 
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For strong dissipation one observes jumps in C/N — meaning that a few particles 
perform many collisions without affecting the global energy too much. In that case one 
expects that the time between collisions drops below the threshold tc and r is set to unity 
according to Eq. (|^. The kinetic energy K(t) is plotted against C/N in Fig. |^(a) for 
simulations with r = 0.60 and different = 10~^° s, 10~^ s, 10~^ s, 10~^ s, and 10"'^ s. For 
t < 50 s the decay of energy is almost independent of tc- For larger times the simulations 
with greater tc loose less energy since more collisions are elastic. From the plot of K{t) 
against C/N in Fig. P(b) one observes for C/N < 100 a power-law behavior K{t) oc [C/NY 
with C ^ 2.5 and for C/N > 100 a much faster decay of energy, with C ~ 7. 

4.2 Cluster growth 

In the following we discuss a simulation with = 79524 particles, L = 500, g = 0.25, and 
tc = 10~^s in more detail. In Fig. p!0|(a) the energy K(t) is plotted against the simulation 
time. The solid line indicates a slope of —2 as one would get in the homogeneous cooling 
regime, see Eq. (pISl). However, the simulation seems to follow, in average, a slope of -1, as 
indicated by the dashed line. Thus the cooling is much slower when it is non-homogeneous. 
Fig. |l3(b) shows that the fraction of particles that have collisions with a collision frequency 
larger than 1/tc is only about 0.1 percent of the total number of particles for t < 50 s. For 
larger times, the number of elastic collisions increases, but never above three percent of 
C/N. 




Figure 10: (a) K{t) as function of time t for a simulation with = 79524, L = 500, 
g = 0.25, and tc = 10~^s. (b) C/N as function of t for the simulation in (a). The circles 
give the mean number of collisions per particle and the squares give the mean number of 
elastic collisions per particle. 



In Fig. |TT] snapshots of the simulation in Fig. |T0| are displayed. With increasing time t 
(and C/N) structures build up in the system and grow in size. In the bright regions in the 
centers of the clusters the collision frequency is largest. 
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t = 0.640 s, C/N = 39 t = 10.24 s, C/N = 183 




Figure 11: ED simulation with = 79524 particles in a system of size L = 500, volume 
fraction g = 0.25, restitution coefficient r = 0.8, and critical collision frequency 1/tc = 10^ 
s~^. The collision frequency is color-coded — red and blue correspond to small and large 
collision frequencies, respectively. 
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4.3 Probability distribution of the collision frequency 

For a more quantitative analysis of the clustering, the probability distribution for particle 
collision frequencies is examined. P{Nc) gives the probability to find a particle that carried 
out Nc collisions per unit time in the last time interval At = t/2. In Fig. |l2|(a) P{Nc) 
is plotted for the homogeneous cooling during short times. The shape of P{Nc) resembles 
a Poisson distribution whose mean and width xj lS.t continuously decrease with time for 
CjN < 39. For C/N = 8, i.e. t = 0.04 s, the dashed line corresponds to Po(C/At,6.5) At 
with the Poisson distribution for the number of collisions C per particle per time At 

P.(C,x) = ^^£t^. (20) 

with mean and width x/^t = 6.5/0.02 s= 325 s^^. This distribution indicates that the 
collisions are uncorrelated events at the beginning of the simulation. 
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Figure 12: Semi-log plot of P{Nc) for dif- 
ferent C/N values from the simulation 
displayed in Figs. ^ and 11 . The system 
is (a) in the homogeneous cooling state, 
(b) in the cluster growth regime, and (c) 
in the long-time state when the largest 
cluster has reached the system size. 
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At longer times, the probability for a large number of collisions can be well fitted by 
an exponential 

i'e(iV„A) = iexp(-^) , (21) 

with the mean collision rate A = 19 in Fig. [T5(a). As soon as clustering occurs, the form 
of the probability of the number of collisions changes dramatically. In Figs. |12|(b) and (c), 
P{Nc) is displayed at different times, i.e. different C/N. The probability for large (small) 
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collision frequencies increases (decreases), and can now be approximated by a power-law of 
the form 

Pp(N„ B) = 2 , (22) 

with the rate B = 8 s~^. All functions P{Nc) are normalized for arbitrary x, A or B, 
but Pp{Nc, B) has a power-law tail so that the first moment diverges. This indicates the 
inelastic collapse, i.e. the divergence of the global collision frequency. However, in the 
simulation presented here, the collision frequency is limited since tc is non-zero. 

In summary, we found that for C/N > 70 the shape of P{Nc) changes from an ex- 
ponential decay for large Nf. to a power-law behavior. This corresponds to the density 
instability when clusters evolve and grow with time, i.e. cooperative motion and interac- 
tion. Interestingly, the shape of the probability distribution is only weakly varying over one 
order of magnitude in time, from C/N ^ 183 to C/N ^ 1404. At these times, particles that 
carry out many collisions coexist with those which carry out only a few. For even larger 
times C/N > 2567 the function P{Nc) changes shape again, as displayed in Fig. [12|(c); the 
probability for large collision frequencies drops again. 



5 Simulating the Flow through Pipes 

In this section we present another collective phenomenon in a totally different geometry. 
Recent observations of approximately V-shaped microcracks in vertically vibrated sand- 
piles p3[ led to the problem of gravity driven vertical motion of sandpiles in 2D pipes ||^ 



In this situation the material is accelerated by the gravitational force and confined 
to the pipe by two side walls. During the fall, cracks develop in the lower part of the pile 
and ascend progressively inside the bulk in both experiment and simulation . For details 
on experiments and simulations see Refs. 

The reasons for a crack to open have been identified as the fluctuations of either the 
wall surface [0 or the random particle motion It was reported that fluctuations lead to 
a momentum wave in the material. A part of the material is decelerated and the material 
from above hits this slower plug and causes a new, possibly stronger momentum wave. 
The increased pressure on the sidewalls may lead to an even stronger deceleration so that 
eventually a crack opens below the plug — and becomes visible only that late. Thus, cracks 
are a rather bad indicator of the dynamical processes occuring inside a granular material, 
since they need some time to open, before they get visible. In simulations one has access 
to quantities that allow deeper insight what is going on in the material. The number of 
collisions per unit time in which a particle participates can be measured and visualized 
easily. 

The system is a box of width L = l/d and initially filled with particles with diameter 
d, situated on a triangular lattice with lattice constant s = 1.01 d. Each particle is assigned a 
random velocity, uniformly distributed in the range —vq < Vi{0) < vq in both, the horizontal 
and vertical directions. This rather regular system is now allowed to reach a steady state, 
i.e. we start the simulation at t = — t^, using r = r^^ = 1 and /i = /x^ = 0. A typical average 
velocity in our simulations is U = y/ < v'^ > = 0.05 ms~^ for t = Os. Due to the rather 
low kinetic energy, the array of particles is still arranged on a triangular lattice, except for 
a few layers at the top which are fiuidized. In a typical simulation, we use L = 20.2 and 
A^ = 1562, so that the array consists of about 80 layers. At time t = Os we remove the 



bottom, switch on dissipation and friction and let the array fall. In Fig. |T3| we present data 
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Figure 13: Snapshots from a simulation with N = 1562, L = 20.2, e = = 0.9, jj, = jd^ = 
0.5, and [3q = = 0.2. Black particles have a collision frequency of Nc = Os~^, white 
particles have > 10^ s~^, and grey particles interpolate between the two extremes. 



from one special simulation with = 1562, L = 20.2, e = = 0.9, jj, = fi^ = 0.5 and 
= Pow = 0.2. Light (dark) particles carried out many (few) collisions in the last interval 
At = 0.001 s. 

The black bar at the right wall marks the particle at which one observes at first a large 
collision frequency N^. at time t = 0.048 s. Already 0.004 s later the neighboring particles 
react and also carry out more and more collisions. This increase in collision frequency leads 
to an increase of pressure that, in return, leads to more friction, a slow-down of the particles 
close to the wall, and to even more collisions with the following particles. The grey region 
indicates an arch-like structure that spans the whole width of the system at t = 0.056 s. 
Again 0.004 s later this region of large Nc and large pressure has almost disappeared, and 
for t > 0.060 s at the position of the initiating particle (at the black bar) a crack is visible. 
Even later, new arches above and below the original one appeared. 

In Fig. |14| we present again the probability P{Nc) to find a particle in the system that 
performed a number of Nc collisions per unit time within the last time interval At = 0.005 s. 
P{Nc) is calculated from the simulation in Fig. |T3| at times t = 0.02 s, 0.03 s, 0.04 s, and 
0.05 s, i.e. before the first arch (or crack) occurs. The probability for the number of collisions 
can be well fitted by the exponential Pe{Nc,A), see Eq. (pi]), with the mean collision rate 
A = 1 s~^. As soon as arches occur, the probability of the number of collisions changes 
shape. In Fig. ^(h), P{Nc) is displayed at times t = 0.055 s, 0.065 s, 0.075 s, and 0.085 s. 
The probability for large (small) collision frequencies increases (decreases), and can now be 
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approximated by the power-law Pp{Nc, B) with the rate B = 14 s^^ from Eq. (^21). Note 
that the curves in Fig. |1^ and Fig. |1^ are both Eq. (^), however, with different B. They 
look different because the horizontal axis is linear in the first but logarithmic in the second 
figure. 
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Figure 14: (a) Semi- log plot of the probability distribution P{Nc) against the number of 
collisions per unit time from the simulation in Fig. |13[ The data correspond to a time before 
cracks and arching occurs, the line gives Pg from Eq. ( |21| ) with A = 1. (b) Log-log plot of 
P{Nc) from the same simulation as in (a) at later times, when cracks and arches exist. The 
solid lines gives Eq. (^) with B = 14. 



6 Summary and Conclusion 

With a rather simple description of a granular material as an ensemble of inelastic spherical 
particles we have investigated interesting effects like the clustering instability and arching. 
A freely cooling granulate builds up a density instability and clusters of particles are formed. 
The structure on small length scales, catched by the event driven simulation, could not be 
reproduced by the stochastic DSMC method. The large scale structures, i.e. the structure 
factor of the clusters, however, was nicely reproduced with the stochastic method. 

DSMC assumes molecular chaos. Thus we examined different systems and found in 
2D a constant probability distribution for the impact parameter, i.e. molecular chaos, when 
the system is elastic or slightly inelastic. Only for large densities, large systems and strong 
dissipation, the molecular chaos assumption becomes invalid. In those cases, we observed 
a shearing motion of grains, a phenomenon that can be found also in large, dilute systems 
after dense clusters have formed. DSMC cannot model dense clusters and also ED has 
problems to handle dense regions with large collision frequencies. Therefore, we introduced 
the advanced TCED method that avoids the inelastic collapse, i.e. dissipation is switched off 
when too much collisions occur per unit time. We found that this variation of the traditional 
ED method involves only a small percentage of the particles and thus does not affect the 
system's behavior as long as the cut-off time tc is small enough. 

Finally, the statistics of the particle collision frequencies was examined for two to- 
tally different boundary conditions. In both cases, freely cooling and pipe flow as well, 
the probability distribution of the collision frequencies shows two types of behavior. In 
the homogeneous, random regime the distribution resembles a Poisson distribution (i.e. an 
exponential for small mean values), indicating that the collisions are uncorrelated events. 
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As soon as cooperative effects like clustering or arching occur, the probabihty distri- 
bution changes from an exponential to a power-law shape. We proposed one functional 
form that approximated the distribution functions for both boundary conditions equally 
well. The described cooperative phenomena are of local nature but leave a fingerprint, 
i.e. a power-law, in the global distribution function of the collision rate. An open issue is 
the theoretical verification and understanding of the shape of this interesting probability 
distribution function. 
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